    #include <math.h>
    #include <time.h>
    #include <mpi.h>

    typedef struct
    {
      float v;

      int x,t;
    } min_Vp;

    #define max_P 8
    #define dim_n 2
    #define max_n 1000 //nodes @ each process

    float frand_Sin(double);
    float distance(float *,float *);

    int main(int argc,char *argv[]) 
    {
      int c1,c2,c3,rank,proc,size,offset=0;
      int cnt=0,stop=max_n*max_P-1;
      min_Vp min;
      float mjn;
      int index[max_P];         //dim_n=8 -> {7,8,9,10,11,12,13,14}
      int edges[max_P*2-2];     //        -> {1,2,3,4,5,6,7,0,0,0,0,0,0,0};
      int point[max_n];         //point in tree
      float recnt[ 2+dim_n];    //most recent point
      float rebuf[(2+dim_n)*max_P]; //recv buffer for P0
      int taken[max_n];         //order of chosen
      int graph[max_n];         //predecessor in tree
      float nodes[dim_n*max_n]; //dim X n
      min_Vp dist[max_n];        //minimum distance trace
      double tic=clock();
      MPI_Comm comm;

      MPI_Init(&argc,&argv);
      MPI_Comm_rank(MPI_COMM_WORLD,&rank);
      MPI_Comm_size(MPI_COMM_WORLD,&size);

      for (c1=0; c1<max_P; c1++) {
        index[c1]=max_P+c1-1;
        if (c1<max_P-1) {
          edges[c1]=c1+1;
          edges[max_P-1+c1]=0;
        }
      }

      MPI_Graph_create(MPI_COMM_WORLD,max_P,index,edges,1,&comm);

      index[0]=1;
      for (c1=0; c1<max_n; c1++) {
        dist[c1].v=1e20;
        point[c1]=0;
      }
      for (c1=0; c1<dim_n*max_n; c1++)
        nodes[c1]=frand_Sin(cos((float)rank+.2));
      if (rank==0) {
        export_init();
        min.x=0; min.t=-1;
        recnt[0]=0; recnt[1]=0;
        memcpy(&recnt[2],nodes,4*dim_n);
      }

      do {
        MPI_Bcast(recnt,2+dim_n,MPI_FLOAT,0,comm); //0->all: send recent point info
        if (rank==(int)floor(cnt/max_n)) //update taken
          taken[cnt%max_n]=recnt[1];
        if (rank==(int)recnt[0]) {       //all: if point in nodes of recent
          point[min.x]=1;                //point in tree now
          graph[min.x]=min.t;            
        }
        mjn=1e20; min.v=1e20; proc=0;
        for (c1=0; c1<max_n; c1++) {     //all: update minimum tracker
          if (!point[c1]) 
         {
            mjn=distance(&nodes[dim_n*c1],&recnt[2]);
            if (dist[c1].v>mjn) {
              dist[c1].v=mjn;
              dist[c1].t=cnt;
            }
            if (dist[c1].v<min.v) {
              min.v=dist[c1].v;
              min.x=c1;
              min.t=dist[c1].t;
              proc=1;
            }
          }
        }
        recnt[0]=min.v;                  //all->0: send local minimum
        recnt[1]=(double)min.x;
        if (proc) memcpy(&recnt[2],&nodes[min.x*dim_n],4*dim_n);
        MPI_Gather(recnt,2+dim_n,MPI_FLOAT,rebuf,2+dim_n,MPI_FLOAT,0,comm);
        if (rank==0) {                   //0: find gobal minimum
          min.v=1e20;
          for (c1=0; c1<max_P; c1++)
            if (rebuf[c1*(2+dim_n)]<min.v) {
              min.v=rebuf[c1*(2+dim_n)];
              proc=c1;
            }
          recnt[0]=(float)proc;
          recnt[1]=rebuf[proc*(2+dim_n)+1]+(float)proc*max_n;
          memcpy(&recnt[2],&rebuf[proc*(2+dim_n)+2],4*dim_n);
        }
        offset+=max_n;
        if (!rank&&!(cnt%1000)) printf(".");
      } while (cnt++!=stop);

      index[0]=max_n; index[1]=dim_n; index[2]=max_n;
      for (c1=0; c1<max_P; c1++) {       //export graph & coordinates
        if (rank==c1) {
          export_flush(1,index,4,taken);
          export_flush(1,index,4,graph);
          export_flush(2,&index[1],6,nodes);
        }
        MPI_Barrier(comm);
      }
      if (!rank) printf("\nelapsed_time =\n   %10.5f\n",(clock()-tic)/CLOCKS_PER_SEC);
      MPI_Finalize();
      return 0;
    }

    float distance(float *vec1,float *vec2) {
      int c1;
      float dist, norm=0;
      for (c1=0;c1<dim_n;c1++) {
        dist=(vec1[c1]-vec2[c1]);
        norm+=dist*dist;
      }
      return sqrt(norm);
    }

    float frand_Sin(double seed) {
      static double drandSin_cnt;
      float t0=sin(++drandSin_cnt+seed)*1e4;
      return t0-floor(t0);
    }


